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, Motivated by recent research of Nikitin et al (J.Phys.D 49,055301(2009)), we examine the effects 

' of interatomic interactions on adatom surface diffusion. By using a mean-field approach in the 

random walk problem, we derive a nonlinear diffusion equation and analyze its solutions. The 
results of our analysis are in good agreement with direct numerical simulations of the corresponding 
discrete model. It is shown that by analyzing a time dependence of adatom concentration profiles 
^ ' one can estimate the type and strength of interatomic interactions. 

+-> \ 

d , I. INTRODUCTION 

> 

Diffusion is ubiquitous in Nature. It determines the behavior and controls the efficiency of many biological and 
technological processes. Examples include the wetting, conductivity of biological membranes, catalysis, growth of 
' crystals, sintering, soldering etc. Surface diffusion is particularly important in nano-technological processes which are 
aimed to obtain objects of submicron sizes where the surface properties are of the same importance as the bulk ones. 
q . Macroscopic description of diffusion is based on Fick's law, which postulates proportionality between the particle 
flux and the concentration gradient. Establishing a link between macroscopic laws of diffusion and microscopic 
non-equilibrium density matrix approach is one of the most challenging and important problems of non-equilibrium 
^vq , statistical mechanics (for reviews see, e.g. [iHs)]). Surface diffusion is essentially many-particle process. Even at very 
low coverage when the interaction between adatoms is negligible, the random walk of an isolated adsorbed particle is 
I a collective motion due its interaction with substrate atoms 0, 0j . The walker moves in a potential landscape which is 
I changed by the walker Q. The walk on a deformable medium when the walker leaves behind a trail and due to slow 
relaxation the trail affects the next walker, is also a collective process @, @]. At finite coverage particles at surfaces, in 
,_H ' addition to the interaction with a substrate, experience lateral interactions of different origin: the attractive van der 
Waals, oscillating electronic exchange and multipole-multipole electrostatic interactions. The dipole-dipole interaction 
which is due to a polar (mainly dipolar) character of adsorption bonds is long-ranged (as r~ 3 ) and generally repulsive 
since all dipole moments are essentially parallel (see review papers, e.g. 

Bowker and King [it]] used Monte-Carlo simulations in order to clarify the effect of lateral interactions of adatoms 
on the shape of evolving concentration profiles in surface diffusion. They showed that the intersection point of the 
diffusion profiles with the initial stepwise profile lies above max /2 in the case of lateral repulsion and below 9 rnax /2 
in the case of attraction (9 m ax is the maximum concentration in the initial step). 



X 



In a quite recent paper ll| an approach based on error function expansion was proposed to fit experimental concen- 
tration profiles. This algorithm provides a high-accuracy fitting and allows extracting the concentration dependence 
of diffusivity from experimental data. The goal of our paper is to model and examine the effects of interatomic interac- 
tions on adatom surface diffusion. Starting with nonlinear random walk equations where the interatomic interactions 
are considered in the mean-field approach, we derive a nonlinear diffusion equation and analyze its solutions. The 
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results of our analysis are in good agreement with direct numerical simulations of the corresponding discrete model. 
It is shown that by analyzing a time dependence of adatom concentration profiles one can estimate the type and 
strength of interatomic interactions. The paper is organized as follows. In Sec.l we present the model. In Sec. 2 we 
study both analytically and numerically interaction effects for the case of low adatom coverage. Sec. 3 is devoted to 
analytical treatment of nonlinear diffusion in the case of non-monotonic concentration dependence of the diffusion 
coefficient. We also compare our results with the results of full scale numerical simulations and results of experimental 
observation. Sec. 4 presents some concluding remarks. 



II. MODEL AND EQUATIONS OF MOTION 



The transport of particles on the surface is described by the set of random walk equations 

j t 9 n {t) = [Wn+p^n On+p{t) (l - 6 H (t)) - W H ^ n+ p9 H {t) (l - 0n+p(t))] M 
p 

where 9 ft is the probability for a particle occupying the n-th binding site on the surface (in the literature on surface 
science this quantity has a meaning of coverage), Wft-^+g gives the rate of the jumps from the binding site ft to a 
neighboring site ft + p (the vector p connects nearest neighbors) . The terms ^1 — 9ft(t)^ in Eqs. (dJ take into account 
the fact that there may be only one adatom at a given site or, in other words, a so-called kinematic interaction. 
The probability with which the particle jumps from site ft to a nearest neighbor ft + p satisfies the detailed balance 
condition, 

where Eft is the binding energy of the particle situated at site ft , f3 — 1/ksT, is the Boltzmann constant and T 
is the temperature of the system. For the transition rates we choose 

W^ ¥ p = wpef" s * (3) 

which corresponds to setting the activation energy for a jump to the initial binding energy. Here wp = vpe~^ Eh (wp — 
W-g) is the jump rate of an isolated particle with standard notations: up is a frequency factor and Eb — E,-% is the 
height of the random walk barrier. Inserting Eq. into Eqs. ([I} we obtain that the random walk of particles on 
the surface is described by a set of equations 

j t 9 H (t) = £ wp [(l - a {t)) 9 H+ p(t) e? E *+? - 9 n {t) (l - 6 R+ p(t)) e^*] . (4) 
p 

In the case when the characteristic size of the particle distribution inhomogeneity is much larger than the lattice 
spacing one can replace 9ft and E^ by the functions 9(f) and E(f) of the continuous variable f and, by expanding the 
functions 0{f + p) and E(f+ p) into a Taylor series, obtain from Eqs. (j4]) that in the continuum approximation the 
transport of particles on the surface is described by the equation of the form 

d t 9 = wv{(v9 + f3 9(l-9)VE^e pE } (5) 

where the notation w — | p 1 wp is used. 

p 

We will study the particle kinetics in the mean field approach when the binding energy E is assumed to be a 
functional of the particle density 9(f, t): E(f) = £(9). In this case Eq. ([5]) takes a form of nonlinear diffusion equation 

d t 0= v{d(6)V6} , (6) 
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where 



D{H) , ( (l + /30(l-0)g)e^ (7) 



is a nonlinear (collective) diffusion coefficient. 



III. SURFACE DIFFUSION AT LOW COVERAGE 

In what follows we restrict ourselves to studying particle distributions spatially homogeneous along the y coordinate: 
9(r,t) = 9(x,t). We assume that initially the particles are step-like distributed 

9(x,0)=9 max H(-x) (8) 

where H(x) is the Heaviside step function. By introducing a centered particle density £(x, t) — ((){x,t) — 

0.5 mox ^ /Ornax, we see that the initial distribution £(x,0) is an odd function of the spatial variable x. It is obvious 
that in the no-interaction case (E = const), when the diffusion equation (jB]) is linear, the antisymmetric character of 
the function is preserved for alH > 0. This means that in the case of noninteracting particles the concentration 

profile for each time moment t passes through the point ^0, ^ . However, interacting diffusing particles exhibit 
quite a different behavior. In 1969, Vedula and one of the authors detected for the first time that concentration profiles 
formed in the process of surface diffusion of thorium on tungsten intersected the initial step-like profile at a point 
lying well above max (see 0, Since then, similar behavior has been found for many electropositive adsorbates 

whose adatoms are known to interact repulsively. A recent example obtained in the case of surface diffusion of Li on 



the Dy-Mo (112) surface was discussed in 11] . It is worth noticing that the above mentioned behavior was observed 



even for rather low coverage: 9 max < 0.3 (see Fig. 2 in ll|). Therefore to explain such a behavior one may assume 



that the binding energy E(r) is linearly dependent on the particle concentration: 

E(r) = E a + [ df V(r- r") 9(f) « E + 9{r) [ df* V(f) (9) 



where Eq is a site energy and V(r — r ) is an interaction parameter which includes all types of lateral interactions. In 
this case the nonlinear diffusion coefficient takes a form 



D(9) = D* [1 +a9(l - 9)j e a& (10) 

where 

D* =we f3E ° = ve-^ E »- E ^ 

is the diffusion coefficient for an isolated particle (or a so-called tracer diffusion coefficient) and the dimensionless 
parameter 



, V = J df V{?) 



characterizes the strength of the lateral interaction. 

The aim of this section is to develop an approach which allows estimating the effects of interparticle interactions 
in the surface diffusion. It is seen from Eqs. (|6]) that the spatio-temporal behavior of the centered particle density 
£(x, t) is governed by the equation 



8r£ = %u+p(.e)), (ii) 
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where t = Do t is a rescaled time and the quantity 

o 

describes the nonlinear properties of the diffusion and vanishes when a — > 0. Taking into account that Eq. ([TTjl with 
the initial condition given by Eq. (j8]) is invariant under gauge transformations r — > A 2 r, x — > A x, 9 — > 9 (A is an 
arbitrary number) one can look for a solution of Eq. (fTl"]) in terms of the Boltzmann variable z = -^p '■ £( x i T ) — C( z ) 
where the function £(z) satisfies the equation 

with the boundary conditions 

C(z)->T~, z^±co. (14) 
Eqs. (fT5|) . (|14l) can be rewritten in the form of the following integral equation 

oo oo 

C(z) = ^ f dzw + {z)-\(l-^- I ' dzwAz))eri{z)- 



= Z 

^ J dzi e£ (erf(z) - erf(z 1 )) ^ p((( 



w ± (z) = (l - erf(z)) [p(c(z)) ± p(t(-zj)] . (15) 

where erf (z) is the error function [l9j . It is seen from Eq. ([T5|) that the concentration profiles {;(x, t) for different time 
moments intersect at the point ^0, C(0)^ with 

oo 

C(0) = # / dzw + (z). (16) 
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In the weak interaction/low coverage limit when a9 max < 1 one can replace the function £(z) in the right-hand-side 
of Eqs. (|15p . (I16p by its expression obtained in the linear case: Coi z ) = \ erf(z) and obtain approximately that under 
the step-like initial condition (JSJ the concentration profiles 9(x,t) for different time moments intersect at the point 
which corresponds to the concentration 



% = (I + C(o)) e v 



({0) = ^-^a9 max ^0.1a9 max . (17) 

47T 

Thus the concentration value 9q at which the concentration profiles intersect changes in the presence of lateral 
interatomic interactions: 9 > 9 max /2 (9q < 9 max /2) when the interaction is repulsive (attractive). We applied 



Eq. dTT)) to analyze the results obtained in llj for the diffusion of Li on the Dy-Mo (112) surface at low coverage 
(dmax ~ 0.33) for which 9q w 0.19 and found out that a9 max w 0.97. It is seen that strictly speaking it is not fully 
legitimate to use our simple analytical perturbation approach (which is valid for a 9 max <C 1) to analyze the results of 



experiments 
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but a qualitative agreement takes place. To validate our analytical results we carried out numerical 
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= M e^"- 1 0jv-i - (l - 0jv_i) e^ £ " 0jv (18) 



simulations of Eqs. (|4|, (|3|) which in the ID-case for a system with N binding sites have the form 

0i = (1 - 0i) e 13 ^ 2 - (l - 2 ) e^ 1 0i 

Jj: n = (1 - 0„) (> f "+i 0„+i + e^- 1 0„_i) - (2 - 0„+i - 0„_i) e^" 0„, (n = 2, ...N - 1) 

d 

where f3 £ n = a0„. Thus, in our model the total number of particles is a conserved quantity. As an initial state we 
used a step-like distribution 

N 

On = Qrnax, for 1 < Tl < — , 

0„ = 0, otherwise . (19) 

We found out that for 6 max = 0.33 the concentration profiles intersect at the point (0, 0.19) ( as it was observed in the 
experiment ll|) for a « (3 ~ 3.5) ( see Fig. [IJ or a0 max « (1 4- 1.5) which is a good agreement with our analytics. 
Thus basing on our approach one can conclude that Li adatoms on the Dy-Mo (112) surface mostly repel each other 
and the intensity of the repulsion is Vq « (3 -j- 3.5) ks T. 



Diffusion of Li adatoms on Dy/Mo(112) was investigated experimentally at T = 600K [ll|, so the estimated 
repulsion energy Vb amounts to 0.16 — 0.18 eV. Let us assess this value in terms of the dipole-dipole interaction. The 
energy of the repulsive interaction between two dipoles having moments p and situated on the surface at a distance 
r is 

Udd = 2 E! « 1-25P 2 [Debyes] 
r 3 r 3 ^Angstroms) 

The dipole moment can be determined from the work function change Atp using the Helmholtz formula for the double 
electric layer: 

| Aip |= innpe , (21) 

where n is the surface concentration of adatoms and e is the electronic charge. For Li on the Mo(112) surface, the p 
value at low coverage was found to be 1.4 Debyes Then, using Eq. (j2T)]) we can find that for two such dipoles 
the interaction energy Udd — 0.16eV can be attained at a distance r w 2.5A, which is close to the distance between 
the nearest adsorption sites (2.73A) within the atomic troughs on Mo(112). This estimation shows that the intensity 
of the lateral interaction deduced from the diffusion data in the way presented above seems physically reasonable. 
Recall, however, that Vq determines a resultant effect experienced by a jumping particle from all its counterparts 
which, in the case of heterodiffusion, are non-uniformly distributed over the surface and provide an additional driving 
force (supplementary to the coverage gradient) that favors a faster diffusion of repulsing particles. 

The adatom interaction effects fade away for the late stage of the evolution when the particle density becomes small 
and the diffusion process transfers into a linear regime ( see Fig. [T]for t = 60000). 



IV. MEAN-SQUARE DEVIATION 

The adatom interactions manifest themselves also in the integral characteristics of kinetics of adatom diffusion. It 
is well known that in the linear regime the variance 

oo 

J dxx 2 9{x,t) 

<- 2 > = ^ (22) 

J dx 6(x, t) 
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FIG. 1: Concentration profile obtained from nonlinear random walk equations (|18[) for a = 3 and initial distribution given by 
a thin dashed curve. The left panel shows the early stage of evolution: wt = 500 (solid curve), wt — 2000 (dashed curve), 
the right panel shows the late stage of evolution: wt — 40000 (dashed), wt = 60000 (solid). The line which corresponds 
to 6 = 9 max /2 is shown as a horizontal dotted line in the left panel and the line which corresponds to the final state of the 
evolution 6 = 6 max /A is shown as a horizontal dotted line in the right panel. 



behaves (in one-dimensional case) as (a; 2 ) = 2t. Therefore it is naturally to introduce a variance rate 

A(t) = (2-±(x 2 )) 2 (23) 

whose time dependence provides a useful information about nonlinear effects in the diffusion process. 
For this quantity we obtain from Eq. P^)l that 



oo oo 




— oo — oo 



Assuming that initially the particles are concentrated in a finite domain in a II-likc form: 

6(x, 0) = 9 max (H{x + I) - H{x - /)) (25) 

where 21 is the size of the initial domain, for small nonlinearities a and low coverage max <C 1 we obtain approximately 
that 



oo oo 



where 




is the solution of the linear diffusion equation with the initial condition (|25p. In the limit of small I we obtain that 

A(r) « " 2 °™ ax I 2 . (28) 

2-7TT 

We checked our analytical considerations by carrying out numerical simulations of Eqs. (fTgf with the initial 
concentration profile given by Eq. (|25|) (see Fig. [2]) for different values of the nonlinearity parameter a. The results 
of these simulations are presented in Fig. [3J The figure shows that the numerically evaluated temporal behavior of 
the rate function A is in a good agreement with our analytical expression given by Eq. ([28)) . Moreover the slopes of 
the curves as it is prescribed by the analytics relate as 0.51 : 0.91 : 2.0 : 3.6 « a\ : a 2 . : a§ : a\ ~ 0.15 2 : 0.2 2 : 0.3 2 : 
0.4 2 . Thus, by measuring the temporal behavior of concentration profiles, it is possible to estimate the strength of 
interatomic interactions. 
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V. CONCENTRATION PROFILES WITH PLATEAU 

In general, the diffusion coefficient is a non-monotonic function of atomic concentration (see e g (sj). There is a 
number of physical reasons which can cause a non-monotonic coverage dependence of the diffusion coefficient. In 
thermodynamic terms, the diffusion flux is proportional to the gradient of chemical potential of adsorbed particles fj, 
which can be written as [l^|.[l5| 

H = fio ~ q(0) + ~ In f t^) (29) 



13 VI - . 

The first term in this equation is the standard chemical potential of the adsorbate, q(6) is the differential heat of 
adsorption and the third term stems from the entropy of mixing of adatoms with the vacant adsorption sites on the 
substrate. (Note that this simplified expression relates only to the first monolayer and does not take into account the 
possibility of formation of the second and next monolayers) . The diffusion coefficient can be represented as a product 

D W = ^5¥*. (30) 

where Dj is a so-called kinetic factor (or jump diffusion coefficient) [f|,Q,[3] and the derivative in the brackets is 
named the thermodynamic factor. In a simplest case when the cross-correlations between the velocities of diffusing 
particles are absent, Dj coincides with the tracer diffusion coefficient D* given by Eq. (10). Inserting Eq. (|29p into 
Eq. ([30]), we get 

D(0) = D j (-^ + Y ^). (31) 

It is seen from Eq. (|3Tj) that any effect which entails a sharp decrease in the heat of adsorption as a function of 
coverage will result in a maximum of the diffusion coefficient in this coverage range For instance, such a situation 
occurs when all energeticall y p rofitable sites at the surface are occupied and adatoms start to fill less favorable sites. 
Actually, Bowker and King [10] found in their Monte Carlo simulations that a well-pronounced maximum in the D{9) 
dependence observed by Butz and Wagner can be explained by existence of two types of lateral interactions: 
repulsive one between the nearest neighbors and attractive between next-nearest neighbors. A similar effect is typical 
for volume diffusion of interstitial atoms in disordered binary alloys having a BCC structure with two nonequivalent 



interstitial positions [21| . In the framework of local equilibrium statistical operator approach |22j it was shown that the 



physical reason for a non-monotonic concentration dependence coefficient is a combined action of lateral interaction 



and adatom density fluctuations [23J. A sharp drop in the heat of adsorption is also observed in the transition 
from filling the first, strongly bound (chemisorbed) monolayer to filling the second, weakly bound (e.g., physisorbed) 
monolayer. In such a case, the spreading of the first monolayer proceeds through diffusion in the mobile uppermost 
(second or next) monolayer (the so called "unrolling carpet" mechanism) [l|. This example shows that a change in 
the heat of adsorption can be accompanied not only by variation of the diffusion parameters (the activation energy 
and prefactor Do in the Arrhenius equation), but also by a change in the atomistic diffusion mechanism itself. 

It is worth noting also that the non-monotonic concentration dependence may be phenomenologically connected 



with a step-like dependence of the heat of adsorption on the coverage ( see a review paper [2J]) . Fig. @] shows the 
diffusion coefficient calculated from Eq. ([7]) by assuming that the on-site adatom energy £{&) (which in most cases is 
proportional to the heat adsorption q(0)) has the form 



£{6) = a(l + tanh/e(0-0 t / l r)J (32) 

where the parameter Q t hr gives the threshold value of the coverage and the parameter k characterizes the sharpness 
of the transition to the new state 25 1. 
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The aim of this section is to consider the diffusion process with a step-like concentration dependent on-site energy 
given by Eq. (|32[) and clarify what kind of new information one can derive by comparing theoretically obtained 
concentration profiles with experimentally obtained ones. 

It is very hard and may be hopeless to solve the equation (j6]) with the diffusion coefficient given by Eqs. and 
(1321) . However, the problem can be solved and some insight into the kinetics can be achieved in the limiting case of 
very sharp energy concentration dependence: k ~ > oo . In this case the diffusion coefficient (JTJ) and (|3"2")l takes the 
form 



D(8) =D*(l + a 5(6 - 9 thr ) +bH{8- 8 thr )^j 



where 



W(l-0tfcr), b = e 2 



1. 



(33) 



(34) 



The nonlinear diffusion equation © with the diffusion coefficient (|3"3")l and the initial condition © has a self-similar 
solution 9(x, r) = Q(z), (z = x/2y / r) which can be presented in the form (see Appendix for a detailed derivation) 



crfc 



'thr ■ 



(x,t) = < 



crfc Izi 

9thr 

rfc ( - 

erfc ^22 

Here the parameters z\ and z 2 are determined by the equations 



I Omax ~ Sthr ) 77 T 

\ / crfc I 22 J 



when x > 2 Z\ yfr 

when — 2 Z2 \fr e a < x < 2 z\ y/r 

when x < —2 z^ e a yfr 



Z-2 



a (I- 6 th r) e Zl erfc(zi) - z\ e 
'thr) e z ? erfc(z 2 ) = 6 thr e z i erfc(zi) 



(35) 

which are obtained from Eqs. (|A9|) taking into account the definition (p4|) . Thus, the concentration profile in the case 
of non-monotonic diffusion coefficient is characterized by existence of a plateau where the concentration of adatoms 
does not depend on the spatial variable x. The length of the plateau £ p = 2^/t(z\ + Z2e a ) increases with time. 
Such a behavior is shown in Fig. [5] The rate with which the length i + p of the plateau increases is determined 
by the nonlinear parameter a and the threshold coverage Qthr- We carried out also numerical simulations of Eqs. 



181) with the the step-like on-site energy £„ 



and as it is seen from Fig. [5] our simple 



y. 1 + tanh (k {0 n - tfcr )j 

model (|33p is in reasonable agreement with numerics. Note that the plateau in the concentration dependence develops 
only on intermediate stage of the evolution. For large enough times the height of the concentration profile becomes 
small, effects of interatomic interactions are negligible and the profile evolutes in accordance with the linear diffusion 
equation (see Fig. [5]). 

It is worth noting that the diffusion of Dy adatoms absorbed by Mo (1 1 2) for the initial coverage 9(x, 0) ~ Q.7H{—x) 
shows a very well pronounced plateau in the concentration profile dependence both on the spatial coordinate for 
different time moments and on the Boltzmann variable ( see Figs. 7 and 8 in ll|). It means that as it is prescribed 
by our analytical considerations the length of the plateau increases as t 1 / 2 . This suggests that our simple analytical 
model may be a useful tool in analyzing an experimentally observed concentration behavior. 



VI. CONCLUSIONS AND DISCUSSION 



In this paper, we have investigated the role of interactions between adatoms in surface diffusion. The problem 
was considered analytically in the mean-field approach. By analyzing discrete nonlinear random walk equations and 
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FIG. 2: Concentration profiles for initial pulse- like distribution (dashed line) and for r = 200 (solid line). The nonlinearity 
parameter a = 0.2. 

corresponding nonlinear diffusion equations with an initial condition in a form of step-like concentration profile, we 
have found that the interactions between adatoms influence significantly the concentration profile development on 
early and intermediate stage of the process. In the case of low coverage the interaction between adatoms makes the 
concentration profile asymmetric: it is shifted to the side of high concentration in the case of repulsive interactions 
and to the side of low concentration for attractive interactions. By calculating the magnitude of the shift one can 
estimate the intensity of lateral interactions between adatoms. At the late stage of kinetics the role of interatomic 
interactions becomes negligible. By studying the nonlinear random walk process which is characterized by a sharp 
maximum in the concentration dependence of the diffusivity, we have found that a well-pronounced plateau develops 
in the concentration profile. The length of the plateau increases in time as t 1 ^ 2 . The height of the plateau Othr 
corresponds to the maximum of the diffusion coefficient which in the frame of our approach corresponds to a sharp 
decrease in the heat of adsorption as a function of coverage. The rate with which the length of the plateau increases 
with time is determined by an amount at which the adsorption heat drops at the threshold coverage Q t hr ■ All above 
mentioned results can be verified experimentally. 
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The nonlinear diffusion equation ([SJ with the diffusion coefficient (1551) and the initial condition ([5]) has a self-similar 
solution 9(x, r) = O(z), (z — x/2y/r) which satisfies the equation 
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FIG. 3: Numerically obtained variance rate A given Eq. (|23|) as a function of the inverse time 1/t for three different values of 
the nonlinearity parameter a: a = 0.15 (dotted line), a = 0.2 (dashed line), a = 0.3 (solid line). 



0.2 0.4 




FIG. 4: Diffusion coefficient for a step-like on-site energy. The parameters used are a = 1, £ = 0.05, 8q = 0.5 



The boundary conditions for Eq. (|Al|) arc 



Q(z) -t 0, for z -cm , 
6(z) 0, for z -> oo . 



(A2) 




FIG. 5: Analytically obtained concentration profile for the diffusion coefficient in the 5-function limit (|32[) with a = 1, 6o = 
0.5, = 1 for three different time moments : r = 0.005 (dotted line), r = 1 (dashed line), r = 5 (solid line). 
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FIG. 6: (color on-line) Numerically obtained concentration profiles in the case of the step function energy dependence given 
by Eq. (p?2")) with 8 = 0.5, a = 1, £ = 0.005 for different time moments: t = (dotted gray line), t = 6000 (dotted line), 
t = 12000 (dashed line), t = 20000 (solid line), t = 40000 (dashed, red, thick line). 

From Eqs. (|A1[) and (1551) we see that the function 



y(e)= / de'z(e') 



(A3) 



where z(O) is an inverse function with respect to Q(z), satisfies the equation 



de 2 



D(ey 



(A4) 



or equivalently, two equations 



d 2 y 1 



dQ 2 y' 
d 2 y _ 1 + 6 
dQ 2 ~ y ' 



augmented by the jump condition 



and the continuity condition 



dy 

de 



_ dy_ 

e=e thr +o dQ 



for 6 < 9 t hr , 
for 9 > 8 thr 



e=e thr -o 2y(0 c ) 



(A5) 



(A6) 



Vtfthr + 0) = y(0 thr - 0) = y{9 c ) 



(A7) 



By integrating Eqs. (|A5|) once, we get 



rfy /oil y( e thr) 



1 d ^J 2 z 2 + ln^ : 



for 6 < 



for 6 > 



Hhr 



'thr 



(A8) 
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where the constants y(9thr), z\ and z 2 satisfy the equations 



%hr = -z ; 7=F eZl erfc(zi) , 



Vl + b (e max - 9thr) = n r e ^ erfc ( z 2) 



y(0 thr ) = —7=* (A9) 
Zl + V 1 + bz 2 



which were obtained from the jump condition (|A6P and the continuity condition (|AT|) . Taking into account the 
definition (|A3[) . we obtain eventually from Eqs. (IA5|) that the concentration profile is determined by the following 
expressions 

erfc(z) 

0(z) = thr — — — - when z>Z\, 
crfc(zi) 



0(z) = 8thr when — yl + b z 2 < z < Z\ , 
erfcf 

0(z) = 9 max - (e max - 9 thr ) \ when z < -z 2 VT+b . (A10) 

V / erfc(z 2 ) 
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